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We present a method that offers perspectives to perform fully antisymmetrized simulations for 
inhomogeneous bulk fermion matter. The technique bears resemblance to classical periodic boundary 
conditions, using localized single-particle states. Such localized states are an ideal tool to discuss 
phenomena where spatial localization plays an important role. The antisymmetrisation is obtained 
introducing a doubly-periodic structure in the many-body fermion wave functions. This results in 
circulant matrices for the evaluation of expectation values, leading to a computationally tractable 
formalism to study fully antisymmetrized bulk fermion matter. We show that the proposed technique 
is able to reproduce essential fermion features in an elegant and computationally advantageous 
manner. 



Bulk fermion systems are ubiquitous in nature. Mostly 
studied as periodic homogeneous structures, they often 
exhibit spatial localizations resulting from impurities, 
random external fields or local minima in the overall po- 
tential [1-3 . The crust of a neutron star is an example 
of such a system. It is build up from protons, neutrons 
and electrons governed by short-range nuclear attraction 
and long-range Coulomb repulsion. As a result of these 
interactions and the subnuclear densities of the neutron 
star's crust, spatial localizations play a crucial role. In 
the lower-density regions of the crust, the crustal mat- 
ter organizes itself in a Coulomb lattice, while in the 
higher-density regions (i.e. the crust-core interface), a 
subtle interplay between those interactions lead to com- 
plex structures dubbed "nuclear pastas." With a prepon- 
derance of low-energy excitation levels, these shapes are 
susceptible to low-energy dynamics stemming from exter- 
nal probes or temperature changes. Molecular dynamics 
techniques are appropriate to study such systems. Of the 
many existing techniques [4]-[8], antisymmetrized molec- 
ular dynamics (AMD) and fermionic molecular dynamics 
(FMD) are the ones that allow a full quantum treatment 
of a fermion system, using an antisymmetrized set of lo- 
calized states. However, a formalism to employ AMD 
or FMD on bulk systems with full antisymmetrization is 
still missing. 

In this paper, we introduce a technique for simulating 
bulk fermion systems that is based on a doubly-periodic 
structure in the many-body wave function. The method 
makes use of localized nonorthogonal single-particle 
states and allows one to address spatial localizations. 
We demonstrate how a complete antisymmetrization 
of localized states under periodic boundary conditions 
can be achieved in a computationally attractive fashion. 
Furthermore, the presented technique is directly appli- 
cable to AMD and FMD. 
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When investigating the properties of bulk matter by 
means of large simulation volumes, the evaluation of ex- 
pectation values becomes cumbersome. Moreover, sur- 
face effects may influence the results when a large frac- 
tion of the constituents lies on the surface of the simu- 
lation volume. A time-honored method is to introduce 
a periodic structure in the simulation. When studying 
bulk matter using a periodic structure, it is crucial to 
make sure that the studied properties of the small but 
infinitely repeated periodic system and the macroscopic 
system which it represents, are the same. As long as the 
correlation volume of the interactions does not exceed the 
simulation volume, the imposed periodicity works fine. 
However, serious problems arise in the presence of long- 
range correlations such as for example those induced by 
Coulomb interactions or by the Pauli exclusion principle. 

When studying fermion systems in a mean-field ap- 
proach, the many-body state is often introduced as an 
antisymmetrized product of single-particle states [9^ . For 
the description of bulk matter, the single-particle states 
of Slater determinants generally fulfill certain boundary 
conditions. Within a periodic structure, the full Hilbert 
space of these single-particle states is spanned by Bloch- 
Floquet states |lQf[12] . These are, however, tedious to 
evaluate and often a subset of the Hilbert space is used. 
This subset is generally referred to as "periodic boundary 
conditions" or the more general "twist-boundary condi- 
tions". The undesired finite-sized shell effects stemming 
from such a restricted set are much reduced when aver- 
aging over the twist angle. This technique is generally 
referred to as "phase-randomization" or "twist- averaged 
boundary conditions" [I3l [14] . 

To describe, however, spatial localizations with delo- 
calized states, a large configuration space is required. To 
circumvent this problem Wannier states can be used as 
single-particle states. Wannier functions span the full 
Hilbert space of the periodic system and represent local- 
ized states [T?,T5lfT6l. 

All previously proposed states fulfill periodic proper- 
ties by definition. In this paper, we propose to represent 
the inhomogeneous bulk system using nonorthogonal lo- 



calized states 0(r), independent of any periodic frame- 
work, and embed them in a periodic-boundary frame- 
work. As the Pauli exclusion principle introduces long- 
range many-body correlations, it poses a major challenge 
to the study of bulk fermion systems with nonorthogo- 
nal localized states. It is indispensable to fully antisym- 
metrize the many-body wave function representing the 
system. 

The properties of a Fermi system are studied by eval- 
uating various A^-body operators. For a many-body 
fermion wave-packet |<l>), denoted as 



many-body fermion Slater determinant is written as 



|$)=i|^l)0---0|^A), 



(1) 



the expectation values of one- and two-body operators 
are calculated as 



^1= ^{(t>p\Bl\4>q)Oq 



(2a) 
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^11 = - ^2 {^P^r\Bll\(t)q(t)s){OqpOsr -OqrOsp). (2b) 

pqrs^l 

Here, the matrix o represents the inverse of the overlap 
matrix n with rvpq = {(j)p\(j)q). Because of the determinant 
structure of |<l>), the expectation values can be written as 
traces using the inverse overlap matrix o. Thereby, the 
matrices n and o play a fundamental role because they 
carry all the information about the fermion statistics of 
the system under study. As their eigenvalues can cover 
many orders of magnitudes, it is paramount to calculate 
these matrices as accurately as possible (i.e. analytically). 
Due to the long-range character of the Pauli correlations, 
this becomes a tedious task for bulk fermion matter. The 
dimension of the matrices usually impedes simulations 
with a large number of particles, however desirable these 
may be for bulk fermion systems. 

Two general approaches exist for creating bulk sys- 
tems : positioning the single-particle states periodically 
or using periodic single-particle wave functions. While 
the first technique leads to infinite matrices, the sec- 
ond results in a limited antisymmetrization. We pro- 
pose a technique which overcomes these problems by in- 
troducing a doubly periodic structure in the description 
of bulk fermionic matter. As can be seen in Fig. [T] 
both the spatial positioning of the simulation volume 
and the single-particle wave packets are made periodic. 
A number of particles are placed in unit cells on a lat- 
tice 05, which tessellate space perfectly. Each cell, con- 
taining A particles, can be identified by a lattice vector 
R = riiai -\- n2a2 -\- nsas with integer rij. A simula- 
tion of order m consists of the cells on the lattice vec- 
tors defined by the condition Uj G {— r?i, . . . ,7ti}, and 
contains {2m + 1)^ identical copies. To eliminate the 
surface effects, the single-particle states are subjected to 
the Born- von Karman boundary conditions [121 [11] with 
periodicity {2m + l)«j. After these manipulations, the 
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where {S}m is the finite set of lattice vectors, T{R) the 
translation operator over R, and 
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The trial state \^rn) can be understood as an infinite 
system with a truncated range of the antisymmetry op- 
erator. Alternatively, \^m) can also be interpreted as a fi- 
nite system with limited periodicity mapped on a toroidal 
structure. In case of the example given in Fig.[l] the orig- 
inal 9 cells (A through F and the central gray cell) are 




FIG. 1. (Color online) Simulation volume (full-bordered cells) 
for a two-dimensional periodic system of first order, created 
by a unit-cell (gray-shaded cell) and eight identical copies (A 
to H) with box length i, symmetrically placed on the lattice 
sites of a two-dimensional square lattice. In addition, the area 
outside the actual simulation volume (dotted area) is gener- 
ated by the particles' wave functions that have a periodicity 
3i in both directions. Therefore, the constituents of boxes 
with identical labels are the same. Because of the double 
periodicity, the constituents in cells on the border of the sim- 
ulation volume not only feel the influence of the particles in 
the neighboring cells but also of all other particles. For in- 
stance, the particles in cell A feel the effect of the particles 
in the central cell and in cells B, C, G and H directly, while 
the periodicity of the wave functions introduces indirect in- 
teractions with the particles in all emanations of cells D, E 
and F. Hence, the double periodicity mimics an infinite sim- 
ulation volume in an effective way. The method is iterative: 
while the two-dimensional simulation volume can be created 
by a single unit-cell, it can also be seen as the periodic ex- 
tension of a one-dimensional periodic simulation volume (the 
blue-shaded row in the figure) created by the unit cell in the 
center. 



mapped onto a torus. The dotted regions are then ob- 
tained by unfolding the torus while keeping the toroidal 
boundary conditions intact [see Eq. Q]. 

Due to the doubly-periodic structure of the trial state 
1^^), the fundamental overlap matrix and its inverse ex- 
hibit a peculiar nested block-circulant structure. This 
circulant structure can be exploited to map the descrip- 
tion of the complete trial state onto a set of equations 
related to only one unit cell. For the overlap matrix of 
a one-dimensional system, denoted as N, the circulant 
structure is given by 



N = 



no n_i • • • ni-m n_m n^ ■ • ■ n2 ni 

ni no n_i • ■ ■ ni_^ n_^ n^ ■ ■ • n2 

: ni no '■ : Tii-m ^-m '■ '■ 

rim-i '■ '■ ■• n_i : '• "• n^ 

n^ Hm-i ■■■ ni no n_i •■■ ni_^ n_^ 

n_m ^m n^-i • • • ni no n_i • ■ • ni_„ 

: n_TO Yim ■ • : ni no ' • : 

n_2 : ■• ■ n^-i : "• "• n_i 

\ n_i n_2 ■ ■ • n_^ n^ T^m-i ■ " " ni no / 



(5) 



The blocks uj^ are Ax A matrices with elements Upq^j^ = 

{4'm,p\T{—kai)\(j)rn,q)' For a two-dimensional periodic 
system, these blocks are (2m + 1)A x (2m + 1)A one- 
dimensional overlap matrices. This is illustrated in 
Fig. [T] This modular structure can be extended to 
higher dimensions, resulting in overlap matrices with 
a nested block-circulant structure. As a result of the 
translation invariance of the system, each A x A block 
of the overlap matrix can be identified unambiguously 
with the lattice vector R that connects the two cells of 
the bra and ket states. The blocks can be evaluated as 

^pq,R = {(t>m,p\T{-R)\(l)m,q) ' 

Genuine antisymmetrization of an infinite system rep- 
resenting bulk matter cannot be achieved in a restricted 
simulation volume. It requires that the equations be eval- 
uated for m ^ oo. Under those circumstances, the wave 
packets |(/>oo,p) have infinite periods and match the origi- 
nal localized single-particle states |0p). The inverse over- 
lap matrix can be obtained through the following scheme 

m- 

Vbz Jbz 

where BZ represents the first Brillouin zone of the lattice 
05 and /c is a vector in this volume. The strength of 
the proposed formalism reveals itself upon evaluating the 
operators in reciprocal space. Normally, the expectation 
values of the infinite fermion system would be calculated 
by means of Eqs. ([2]), resulting in infinite sums over the 
block structure : 



Re^pq=l 



^p\Bif{R)\^q)OR^ 



These sums, however, translate into integrals over the 
first Brillouin zone which are computationally straight- 
forward to evaluate. The expectation value of a one-body 
operator per unit-cell volume can then be computed as 



^p,i 



T7- [ f2 1^i,pq{k)0,p{k)d^k, (6a) 
Vbz Jbz p^^ 



^i,pM= E((/>p|^/T(i^)|(/>,)e^'^•^. 



(6b) 



Re^B 



Here we assume that operator Bj commutes with the 
translation operator T{R). The two-body operator has 
an analogous structure. 

From a computational perspective, a proper choice of 
the localized single-particle states helps speeding-up the 
computations. Gaussian wave packets come with the 
interesting feature that most matrix-elements can be 
evaluated analytically. Furthermore, the computational 
order for the calculation of expectation values is the 
same for bulk fermion systems as for finite-sized systems 
{N^ for two-body interactions), even though the former 
represent a system with an infinite number of particles. 
On the other hand, the equations for the bulk systems 
involve an extra integration over the Brillouin zone of 
the imposed lattice, increasing the computational effort. 
However, it has been shown that integrals over the 
Brillouin zone can be performed using only a limited 
number of specific k points [19l [20] . 
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In the following, we will show that the proposed tech- 
nique reproduces the features intrinsic to the fermionic 
behavior of the system. The single-particle states are 
Gaussian wave packets of the form {x\ab) = exp{ — (a? — 
6)^ /(2a)} where the complex vector b represents the 
mean position in phase space and a is a complex pa- 
rameter connected with the width of the wave packet. In 
Ref. [7 it was shown that for a hundred periodically po- 
sitioned wave packets, Eqs. (pi reproduce the momentum 
and spatial densities intrinsic to one-dimensional Fermi 
systems. In Fig. [2J we show that our method also repro- 
duces this result by placing a single particle in a unit cell 
of size i. It is clearly visible that with increasing y/a/£ 
and thus with increasing overlap of neighboring states 
and hence increasing influence of the antisymmetrization, 
the fermionic behavior becomes apparent. For ^/a/i = 1, 
the momentum distribution does not show a sharp cut- 
off near the Fermi momentum because the periodic set 
of Gaussian wave packets does not constitute a complete 
set. For ^/a/i -^ oo , the antisymmetry operator projects 
the Gaussian wave packets on plane waves and a sharp 
cutoff is reached. Even though this one-dimensional ex- 
ample clearly hints at free-fermion behavior, it has to be 
assured that fermion-like properties of the simulations 
are not an artifact of the symmetry imposed by the lat- 
tice. The latter is, however, the case and it can be shown 
that the momentum distribution evolves toward the first 
Brillouin zone of the lattice 18 . This is demonstrated in 
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FIG. 2. (Color online) Spatial density (left panels) and mo- 
mentum density (right panels) of equidistant Gaussian wave 
packets with one-dimensional periodicity. One Gaussian with 
variance a per unit cell of size i. For small overlap (top : 
\/^/i = 0.2) the system behaves like one of distinguishable 
particles while for large overlap (bottom : y/a/i = 1) the 
fermionic behavior becomes manifest. For comparison, the 
spatial distribution of a uniform Fermi gas {y/a/i -^ oo) 
is presented (dashed line). The momentum distribution of 
distinguishable particles is also presented (dotted line). The 
value of kp is given by tt/I. 



Fig. |3) where a single wave packet is placed on a hexag- 
onal lattice. With increasing overlap, the momentum 
distribution reflects the first Brillouin zone of the hexag- 
onal lattice. This effect can be expected as the width of 
a single wave packet becomes much larger than the size 
of its unit cell, and the system behaves as plane waves in 
a periodic structure. 

As we introduced periodic boundary conditions to 
eliminate surface effects, a successful investigation of bulk 
matter obviously requires that influences of the geome- 
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FIG. 3. (Color online) The momentum distribution of a single 
fermion under hexagonal boundary conditions. With increas- 
ing overlap the shape of the momentum distribution evolves 
toward that of the first Brillouin zone. The value of ki is given 

by 7v/i. 
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FIG. 4. (Color online) Comparison of the scaled spatial den- 
sity distributions pxi'^ /N of A^ = 25 Gaussian wave packets 
under different conditions. The wave packets are randomly 
placed in the unit cell of a square lattice with unit length i 
and have a variance given by y/a — 0.2^. The mean momen- 
tum of the individual wave packets is set to zero. The densi- 
ties are normalized and coordinates are expressed in units of 
i. The upper-two panels depict distinguishable particles. The 
lower- two panels represent the antisymmetrised case. The left 
and right panels show the effect of a simulation without and 
with periodicity introduced, respectively, as described in the 
text. The black circles show the positions of the centroids of 
the wave packets. Note the difference in color scaling. 



try of the chosen boundary conditions are negligible. In 
Figs. |4] and [5] we investigate the effect of periodicity and 
antisymmetry on the spatial and momentum distribu- 
tions for a situation where 25 single-particle states are 
randomly distributed in a square unit cell. The figures 
demonstrate that the spatial and momentum distribu- 
tion of a periodic antisymmetrised system (lower right 
panels) reflect that of a free fermion system. The per- 
turbations on the Fermi-sphere in momentum space, as 
well as those on the uniform distribution in coordinate 
space can be understood as a result of the local clustering. 
Hence, a simulation with the proposed periodic structure 
reproduces the intrinsic bulk fermion behavior, indepen- 
dent of the imposed periodic structure. When compar- 
ing the density distributions of the distinguishable par- 
ticles with the indistinguishable fermions, it is evident 
that the distributions are quite different. This clearly 
shows that antisymmetry introduces Fermi motion, re- 
distributing the densities. The effect of Pauli repulsion 
becomes visible when comparing the maxima of the den- 
sity distribution with the location of the centroids in the 
non-periodic antisymmetrised case (lower right panel of 
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FIG. 5. (Color online) Comparison of the scaled momentum 
density distributions pkhJ/N oi N — 2^ Gaussian wave pack- 
ets under different conditions. The wave packets are randomly 
placed in the unit cell of a square lattice with unit length I 
and have a variance given by ^/a — 0.2^. The mean momen- 
tum of the individual wave packets is set to zero. The den- 
sities are normalized and coordinates are expressed in units 
of ki — V^tt/^. The upper- two panels depict distinguish- 
able particles. The lower- two panels represent the antisym- 
metrised case. The left and right panels show the effect of 
a simulation without and with periodicity, respectively. The 
black circle represents the Fermi "sphere" of a system of free 
particles with the same mean density, kp — \/ ^.ixN jll. Note 
the difference in color scaling. 



Fig. pi). When enabling periodicity, the effect is small for 
the distinguishable particles, but crucial for a good de- 



scription of the bulk fermion system. This clearly shows 
that antisymmetry, a long-range many-body correlation, 
should not be ignored or truncated for the bulk descrip- 
tion. 

In summary, we propose a technique that allows one 
to simulate fully antisymmetrized inhomogeneous bulk 
fermion matter by means of localized single-particle 
states. A doubly-periodic structure imposed on the 
many-body fermion wave function replaces the classi- 
cal idea of periodic boundary conditions, treating each 
single-particle state as an individual state, but with a 
truncated antisymmetry. This leads to a circulant struc- 
ture in the matrix formalism for the evaluation of ex- 
pectation values. This circulant structure is exploited to 
reduce the numerical cost of the simulation. An impor- 
tant merit of the proposed technique is that it becomes 
possible to perform antisymmetrized calculations in the 
limit of infinite volumes. The structure of the resulting 
equations resembles those of a finite fermion system, but 
requires an extra integration over the first Brillouin zone, 
reflecting the periodic boundary conditions. Thus, the 
truncation in the antisymmetry vanishes, yielding a for- 
malism that allows one to compute expectation values of 
various operators for bulk fermion matter in a computa- 
tionally feasible fashion. Although the equations only ad- 
dress a finite number of particles in a unit cell, they keep 
track of the Fermi statistics of the infinite system. We 
evaluate the validity of this description with a study of 
the spatial and momentum distribution of various lattice 
systems [18 and show that our simulations convincingly 
reproduce free fermion behavior. This shows that full 
antisymmetrization can be imposed on the bulk system, 
a feature that has been considered hard to reach. Fur- 
thermore, the use of localized states makes the technique 
suitable to study inhomogeneous bulk fermion matter. 
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